# This file is part of Sonedyan.
#
# Sonedyan is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation;
# either version 3 of the License, or (at your option) any
# later version.
#
# Sonedyan is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public.
# If not, see <http://www.gnu.org/licenses/>.
#
# Copyright (C) 2009-2012 Jimmy Dubuisson <jimmy.dubuisson@gmail.com>

source("misc.R")

# get all cycles (up to specified length) of g according to Johnson's algorithm
#
# Article describing Johnson's algorithm (1975)
# http://dutta.csc.ncsu.edu/csc791_spring07/wrap/circuits_johnson.pdf
# 
# Python implementation (using networkx library):
# http://networkx.lanl.gov/_modules/networkx/algorithms/cycles.html
#
# @return hash (cycle length -> vector of cycles)
get.johnson.cycles <- function(g, maxdepth)
{
	blocked <- hash()
	B <- hash()
	path <- cstack()
	result <- hash()
	
	# unblock specified vertex and all its ancestors
	unblock <- function(vid)
	{
	  #loginfo(paste("Unblocking node:", vid))
	  #bl <- try(blocked[[as.character(vid)]], TRUE)
	  #if (bl == TRUE)
	  if (blocked[[vid]])
	  {
	    blocked[vid] <- FALSE
	    #bs <- try(B[[as.character(vid)]], TRUE)
	    bs <- B[[vid]]
	    while (length(bs) > 0)
	    {
	      nvid <- bs[[1]]
	      bs <- bs[-1]
	      B[vid] <- bs
	      unblock(nvid)
	    }
	  }
	}
	
	# get graph simple cycles
	# component: subgraph of the initial graph
	circuit <- function(thisnode, startnode, component, depth)
	{
	  closed <- FALSE
	  # get node ids
	  sid = V(component)[startnode]$id
	  vid = V(component)[thisnode]$id
	
	  if (depth > maxdepth)
	  {
	    # loginfo("Max depth was reached")
	    return(TRUE)
	  }
	
	  # loginfo(paste("Pushing and blocking node:", vid, "(startnode is:", sid, ")"))
	  path$push(vid)
	  blocked[vid] = TRUE
	
	  for (nextnode in neighbors(component, thisnode, mode = "out"))
	  {
	    nid = V(component)[nextnode]$id
	    # loginfo(paste("Looking for successor:", nid, "(thisnode is:", vid, "| startnode is:", sid, ")"))
	    if (nextnode == startnode)
	    {
			len <- length(path$.data)
			sc <- paste(path$.data, collapse = '.')
			# loginfo(paste("A simple cycle was found:", sc))
			if (is.null(result[[as.character(len)]]))
		  		result[len] <- sc
			else
		  		result[len] <- append(result[[as.character(len)]], sc)
			closed <- TRUE
	    }
	    else if (!blocked[[nid]])
	    {
	      # loginfo("Recursive call...")
	      if (circuit(nextnode, startnode, component, depth + 1))
	      {
			closed <- TRUE
	      }
	    }
	  } 
	
	  if (closed)
	  {
	    unblock(vid)
	  }
	  else
	  {
	    for (nextnode in neighbors(component, thisnode, mode = "out"))
	    {
	      nid = V(component)[nextnode]$id
	      if (!vid %in% B[[nid]])
	      {
				#loginfo(paste("Adding thisnode:", vid, "to the list of parents of nextnode:", nid))
				B[nid] <- append(B[[nid]], vid)
	      }
	    }
	  }
	
	  path$pop()
	  closed
	}
		
	# V(g) vids are sorted by default
	for (s in V(g))
	{
	  loginfo("--------------------")
	  loginfo(paste("s=", s))
	  
	  # get subgraph induced by nodes above s
	  sub <- induced.subgraph(g, which(V(g) >= s))
	  loginfo(paste("subgraph size:", length(V(sub))))
	  
	  # extract SCCs from subgraph
	  sccs = clusters(sub, mode = "strong")
	  loginfo(paste("# SCCs:", sccs$no))
	  
	  # get SCC containing the smallest vertex id
	  minSccId = sccs$membership[1]
	  loginfo(paste("min SCC id:", minSccId))
	  
	  # get the corresponding SCC from subgraph
	  component = induced.subgraph(sub, which(sccs$membership == minSccId))
	  loginfo(paste("component size:", length(V(component))))
	
	  if (length(V(component)) > 1)
	  {
		# nodes get renumbered in igraph so that the least vertex is always 0
	    #startnode = min(V(component))
	    startnode = 0
	
		# initialize component vertices
	    for (v in V(component))
	    {
	      vid = V(component)[v]$id
	      blocked[vid] <- FALSE
	      B[vid] <- c()
	    }
	
	    circuit(startnode, startnode, component, 1)
	  }
	}
	
	result
}

##########
## custom implementation
##########

# get cycles of length 'l'
get.cycles <- function(g, l)
{
	counter <- 0
	size <- length(V(g))
	
	cycles <- c()
	cnames <- c()
	eids <- c()
	vids <- c()
	# vertex and edge cycle multiplicity
	vm <- c()
	em <- c()
	
	# while graph has more than 1 vertex
	while(length(V(g) > 1))
	{
		# get cycles of length l passing through least vertex (always 1 in igraph)
		cs <- get.vertex.cycles(g, 1, l)
		
		if (length(cs) > 0)
		{
			for (c in cs)
			{
				# we necessarily have a 'base cycle' as we started from the least vertex
				cycle <- as.vector(sapply(strsplit(c, "\\."), as.numeric))
				
				# add the # of removed nodes to the cycle array (vertices get renumbered at each loop)
				cycle.id <- paste(cycle + counter, collapse = ".")
				cycle.name <- tolower(paste(V(g)[cycle]$id, collapse = "|"))
				
				cycles <- append(cycles, cycle.id)
				cnames <- append(cnames, cycle.name)
				
				# a new cycle was added, update vertex and edge multiplicity vectors
				# complete vector of nodes with starting node
				cycle <- append(cycle, cycle[1])
				
				for (i in c(1:(length(cycle) - 1)))
				{
					vid1 <- V(g)[cycle[i]]$id
					vid2 <- V(g)[cycle[i + 1]]$id
					edge <- paste(vid1, "|", vid2, sep = "")
					
					eids <- append(eids, edge)
					vids <- append(vids, vid1)
					vids <- append(vids, vid2)
					
					if (is.null(em) || is.na(em[edge])){ em[edge] <- 1 }
					else { em[edge] <- em[edge] + 1 }
					
					if (is.null(vm) || is.na(vm[vid1])){ vm[vid1] <- 1 }
					else { vm[vid1] <- vm[vid1] + 1 }
				}
			}
		}
		
		counter <- counter + 1
		g <- delete.vertices(g, 1)
		loginfo(paste("Step ", counter, "/", size))
	}
	
	vids <- unique(vids)
	eids <- unique(eids)
	
	loginfo(paste("# vertices ", l, ":", length(vids)))
	loginfo(paste("# edges ", l, ":", length(eids)))
	loginfo(paste("# cycles ", l, ":", length(cycles)))
	
	# vids: vector of vertex names
	# eids: vector of edge names
	# cycles: vector of cycles ids (concatenated set of vertex ids)
	# cnames: vector of cycle names (concatenated set of vertex names)
	# vcm: named vector of k-cycle multiplicities (vertex name -> multiplicity)
	list(v = vids, e = eids, c = cycles, cn = cnames, vcm = vm, ecm = em)
}

# get all shortest cycles (All Nodes Shortest Cycles (ANSC) problem)
get.all.shortest.cycles <- function(g)
{
	counter <- 0
	size <- length(V(g))
	
	cycles <- c()
	cnames <- c()
	eids <- c()
	vids <- c()
	
	# while graph has more than 1 vertex
	while(length(V(g) > 1))
	{
		# get shortest cycles passing through least vertex (always 1 in igraph)
		cs <- get.vertex.shortest.cycles(g, 1)
		
		if (length(cs) > 0)
		{
			for (c in c(1:length(cs)))
			{
				# we necessarily have a 'base cycle' as we started from the least vertex
				cycle <- cs[[c]]
				# add the # of removed nodes to the cycle array (vertices get renumbered at each loop)
				cycle.id <- paste(cycle + counter, collapse = ".")
				cycle.name <- tolower(paste(V(g)[cycle]$id, collapse = "|"))
				
				cycles <- append(cycles, cycle.id)
				cnames <- append(cnames, cycle.name)
				
				# if a new cycle was added, update vertex and edge multiplicity vectors
				# complete vector of nodes with starting node
				cycle <- append(cycle, cycle[1])
				
				for (i in c(1:(length(cycle) - 1)))
				{
					vid1 <- V(g)[cycle[i]]$id
					vid2 <- V(g)[cycle[i + 1]]$id
					edge <- paste(vid1, "|", vid2, sep = "")
					
					eids <- append(eids, edge)
					vids <- append(vids, vid1)
					vids <- append(vids, vid2)
				}
			}
		}
		
		counter <- counter + 1
		g <- delete.vertices(g, 1)
		loginfo(paste("Step ", counter, "/", size))
	}
	
	vids <- unique(vids)
	eids <- unique(eids)
	
	loginfo(paste("# vertices :", length(vids)))
	loginfo(paste("# edges :", length(eids)))
	loginfo(paste("# cycles :", length(cycles)))
	
	# vids: vector of vertex names
	# eids: vector of edge names
	# cycles: vector of cycles ids (concatenated set of vertex ids)
	# cnames: vector of cycle names (concatenated set of vertex names)
	list(v = vids, e = eids, c = cycles, cn = cnames)
}

# get all shortest cycles of the specified length (All Nodes Shortest Cycles (ANSC) problem)
get.shortest.cycles <- function(g, l)
{
	counter <- 0
	size <- length(V(g))
	
	cycles <- c()
	cnames <- c()
	eids <- c()
	vids <- c()
	# vertex and edge cycle multiplicity
	vm <- c()
	em <- c()
	
	# while graph has more than 1 vertex
	while(length(V(g) > 1))
	{
		# get shortest cycles passing through least vertex (always 1 in igraph)
		cs <- get.vertex.shortest.cycles(g, 1)
		
		if (length(cs) > 0)
		{
			for (c in c(1:length(cs)))
			{
				# we necessarily have a 'base cycle' as we started from the least vertex
				cycle <- cs[[c]]
				
				if (length(cycle) == l)
				{
					# add the # of removed nodes to the cycle array (vertices get renumbered at each loop)
					cycle.id <- paste(cycle + counter, collapse = ".")
					cycle.name <- tolower(paste(V(g)[cycle]$id, collapse = "|"))
					
					cycles <- append(cycles, cycle.id)
					cnames <- append(cnames, cycle.name)
					
					# if a new cycle was added, update vertex and edge multiplicity vectors
					# complete vector of nodes with starting node
					cycle <- append(cycle, cycle[1])
					
					for (i in c(1:(length(cycle) - 1)))
					{
						vid1 <- V(g)[cycle[i]]$id
						vid2 <- V(g)[cycle[i + 1]]$id
						edge <- paste(vid1, "|", vid2, sep = "")
						
						eids <- append(eids, edge)
						vids <- append(vids, vid1)
						vids <- append(vids, vid2)
						
						if (is.null(em) || is.na(em[edge])){ em[edge] <- 1 }
						else { em[edge] <- em[edge] + 1 }
						
						if (is.null(vm) || is.na(vm[vid1])){ vm[vid1] <- 1 }
						else { vm[vid1] <- vm[vid1] + 1 }
					}
				}
			}
		}
		
		counter <- counter + 1
		g <- delete.vertices(g, 1)
		loginfo(paste("Step ", counter, "/", size))
	}
	
	vids <- unique(vids)
	eids <- unique(eids)
	
	loginfo(paste("# vertices :", length(vids)))
	loginfo(paste("# edges :", length(eids)))
	loginfo(paste("# cycles :", length(cycles)))
	
	# vids: vector of vertex names
	# eids: vector of edge names
	# cycles: vector of cycles ids (concatenated set of vertex ids)
	# cnames: vector of cycle names (concatenated set of vertex names)
	# vcm: named vector (vertex name -> multiplicity)
	list(v = vids, e = eids, c = cycles, cn = cnames, vcm = vm, ecm = em)
}

# get set of shortest cycles passing through vertex 'v'
#
# g: graph
# v: starting node
get.vertex.shortest.cycles <- function(g, v)
{
	# get the set of nodes pointing to v
	nv <- neighbors(g, v, mode = "in")
	# loginfo(paste("Set of parent nodes of v=", v, ":", paste(nv, collapse = ",")))
	# get all shortest path from node v to node nv
	get.all.shortest.paths(g, v, nv, mode = "out")
}

# get set of cycle path of length 'l' passing through vertex 'v'
#
# g: graph
# sn: starting node
# l: length of the cycle
get.vertex.cycles <- function(g, sn, l)
{
	queues <- hash()
	qkeys <- cstack()
	hpaths <- hash()
	klength <- hash()
	
	paths <- c()
	
	# push starting node
	q <- cstack()
	q$push(sn)
	
	# get the set of children of starting node
	nv <- neighbors(g, sn, mode = "out")
	
	qkeys$push(as.character(sn))
	hpaths[as.character(sn)] <- c(sn)
	klength[sn] <- 0
	nq <- cstack()
	nq$push(nv)
	queues[sn] <- nq
	
	while (length(qkeys$.data) > 0)
	{
		qkey <- qkeys$pop()
		hpath <- hpaths[[qkey]]
		queue <- queues[[qkey]]
		cl <- klength[[qkey]]
				
		while (length(queue$.data) > 0)
		{
			cn <- queue$pop()
			nkey <- paste(c(qkey, cn), collapse = ".")
			npath <- append(hpath, cn)
			nl <- cl + 1
			
			# if the path has the desired length
			if (nl == l)
			{
				# if current node is starting node, we have a cycle
				if (cn == sn)
				{
					#loginfo(paste("Cycle found: ", qkey))
					paths <- append(paths, qkey)
				}	
			}
			else
			{
				# if current node is not one of the path nodes, let's continue
				if (!(cn %in% hpath))
				{
					# get the set of children of current node
					nv <- neighbors(g, cn, mode = "out")

					qkeys$push(nkey)
					hpaths[nkey] <- npath
					klength[nkey] <- nl
					nq <- cstack()
					nq$push(nv)
					queues[nkey] <- nq
				}	
			}
			
		}	
	}
	
	paths
}

# get graph cycles up to the 
# NB: length should be >= 2
print.graph.cycles <- function(name, core, l)
{
	# get set of unlabeled edge indexes
	get.unlabeled.edge.indexes <- function(sub, eids)
	{
		reids <- setdiff(E(sub)$id, eids)
		E(sub)[E(sub)$id %in% reids]
	}
	
	# print cycle labels
	print.cycles <- function(g, cids)
	{	
		for (i in cids)
		{
			for (j in strsplit(i, "\\."))
			{
				print(V(g)[as.numeric(j)]$id)
			}
		}
	}
	
	print(paste("Getting cycles of", name, "core..."))
	
	# 2 cycles
	cycles2 <- get.cycles(core, 2)
	
	# 3 cycles
	if (l >= 3)
		cycles3 <- get.cycles(core, 3)
	
	# 4 cycles
	if (l >= 4)
		cycles4 <- get.cycles(core, 4)
	
	# 5 cycles
	if (l >= 5)
		cycles5 <- get.cycles(core, 5)
	
	# 6 cycles
	if (l >= 6)
		cycles6 <- get.cycles(core, 6)
	
	# cycle ids
	cids2 <- cycles2$c
	
	print("Printing 2-length cycles...")	
	print.cycles(core, cids2)
	
	if (l >= 3)
	{
		cids3 <- cycles3$c
		print("Printing 3-length cycles...")	
		print.cycles(core, cids3)
	}
	
	if (l >= 4)
	{
		cids4 <- cycles4$c
		print("Printing 4-length cycles...")	
		print.cycles(core, cids4)
	}
	
	if (l >= 5)
	{
		cids5 <- cycles5$c
		print("Printing 5-length cycles...")	
		print.cycles(core, cids5)
	}
	
	if (l >= 6)
	{
		cids6 <- cycles6$c
		print("Printing 6-length cycles...")	
		print.cycles(core, cids6)
	}
	
	# vertex sets
	vids2 <- cycles2$v
	
	if (l >= 3)
		vids3 <- cycles3$v
	if (l >= 4)
		vids4 <- cycles4$v
	if (l >= 5)
		vids5 <- cycles5$v
	if (l >= 6)
		vids6 <- cycles6$v
	if (l >= 3)
		vids23 <- union(vids2, vids3)
	if (l >= 4)
		vids234 <- union(vids23, vids4)
	if (l >= 5)
		vids2345 <- union(vids234, vids5)
	if (l >= 6)
		vids23456 <- union(vids2345, vids6)
	
	print(paste("# vertices2:", length(vids2)))
	
	if (l >= 3)
		print(paste("# vertices3:", length(vids3)))
	if (l >= 4)
		print(paste("# vertices4:", length(vids4)))
	if (l >= 5)
		print(paste("# vertices5:", length(vids5)))
	if (l >= 6)
		print(paste("# vertices6:", length(vids6)))
	if (l >= 3)
		print(paste("# vertices23:", length(vids23)))
	if (l >= 4)
		print(paste("# vertices234:", length(vids234)))
	if (l >= 5)
		print(paste("# vertices2345:", length(vids2345)))
	if (l >= 6)
		print(paste("# vertices23456:", length(vids23456)))
	
	sub2 <- get.subgraph(core, vids2)
	
	if (l >= 3)
		sub3 <- get.subgraph(core, vids3)
	if (l >= 4)
		sub4 <- get.subgraph(core, vids4)
	if (l >= 5)
		sub5 <- get.subgraph(core, vids5)
	if (l >= 6)
		sub6 <- get.subgraph(core, vids6)
	if (l >= 3)
		sub23 <- get.subgraph(core, vids23)
	if (l >= 4)
		sub234 <- get.subgraph(core, vids234)
	if (l >= 5)
		sub2345 <- get.subgraph(core, vids2345)
	if (l >= 6)
		sub23456 <- get.subgraph(core, vids23456)
	
	print(paste("# subgraph2 vertices / edges:", length(V(sub2)), "/", length(E(sub2))))
	
	if (l >= 3)
		print(paste("# subgraph3 vertices / edges:", length(V(sub3)), "/", length(E(sub3))))
	if (l >= 4)
		print(paste("# subgraph4 vertices / edges:", length(V(sub4)), "/", length(E(sub4))))
	if (l >= 5)
		print(paste("# subgraph5 vertices / edges:", length(V(sub5)), "/", length(E(sub5))))
	if (l >= 6)
		print(paste("# subgraph6 vertices / edges:", length(V(sub6)), "/", length(E(sub6))))
	if (l >= 3)
		print(paste("# subgraph23 vertices / edges:", length(V(sub23)), "/", length(E(sub23))))
	if (l >= 4)
		print(paste("# subgraph234 vertices / edges:", length(V(sub234)), "/", length(E(sub234))))
	if (l >= 5)
		print(paste("# subgraph2345 vertices / edges:", length(V(sub2345)), "/", length(E(sub2345))))
	if (l >= 6)
		print(paste("# subgraph23456 vertices / edges:", length(V(sub23456)), "/", length(E(sub23456))))
	
	# edge sets
	eids2 <- cycles2$e
	
	if (l >= 3)
		eids3 <- cycles3$e
	if (l >= 4)
		eids4 <- cycles4$e
	if (l >= 5)
		eids5 <- cycles5$e
	if (l >= 6)
		eids6 <- cycles6$e
	if (l >= 3)
		eids23 <- union(eids2, eids3)
	if (l >= 4)
		eids234 <- union(eids23, eids4)
	if (l >= 5)
		eids2345 <- union(eids234, eids5)
	if (l >= 6)
		eids23456 <- union(eids2345, eids6)
	
	reids2 <- get.unlabeled.edge.indexes(sub2, eids2)
	
	if (l >= 3)
		reids3 <- get.unlabeled.edge.indexes(sub3, eids3)
	if (l >= 4)
		reids4 <- get.unlabeled.edge.indexes(sub4, eids4)
	if (l >= 5)
		reids5 <- get.unlabeled.edge.indexes(sub5, eids5)
	if (l >= 6)
		reids6 <- get.unlabeled.edge.indexes(sub6, eids6)
	if (l >= 3)
		reids23 <- get.unlabeled.edge.indexes(sub23, eids23)
	if (l >= 4)
		reids234 <- get.unlabeled.edge.indexes(sub234, eids234)
	if (l >= 5)
		reids2345 <- get.unlabeled.edge.indexes(sub2345, eids2345)
	if (l >= 6)
		reids23456 <- get.unlabeled.edge.indexes(sub23456, eids23456)
	
	print(paste("# reids2 to be removed:", length(reids2)))
	
	if (l >= 3)
		print(paste("# edges3 to be removed:", length(reids3)))
	if (l >= 4)
		print(paste("# reids4 to be removed:", length(reids4)))
	if (l >= 5)
		print(paste("# reids5 to be removed:", length(reids5)))
	if (l >= 6)
		print(paste("# reids6 to be removed:", length(reids6)))
	if (l >= 3)
		print(paste("# reids23 to be removed:", length(reids23)))
	if (l >= 4)
		print(paste("# reids234 to be removed:", length(reids234)))
	if (l >= 5)
		print(paste("# reids2345 to be removed:", length(reids2345)))
	if (l >= 6)
		print(paste("# reids23456 to be removed:", length(reids23456)))
	
	# removing unlabeled edges
	subb2 <- delete.edges(sub2, reids2)
	
	if (l >= 3)
		subb3 <- delete.edges(sub3, reids3)
	if (l >= 4)
		subb4 <- delete.edges(sub4, reids4)
	if (l >= 5)
		subb5 <- delete.edges(sub5, reids5)
	if (l >= 6)
		subb6 <- delete.edges(sub6, reids6)
	if (l >= 3)
		subb23 <- delete.edges(sub23, reids23)
	if (l >= 4)
		subb234 <- delete.edges(sub234, reids234)
	if (l >= 5)
		subb2345 <- delete.edges(sub2345, reids2345)
	if (l >= 6)
		subb23456 <- delete.edges(sub23456, reids23456)
	
	print(paste("# new subgraph2 vertices / edges:", length(V(subb2)), "/", length(E(subb2))))
	
	if (l >= 3)
		print(paste("# new subgraph3 vertices / edges:", length(V(subb3)), "/", length(E(subb3))))
	if (l >= 4)
		print(paste("# new subgraph4 vertices / edges:", length(V(subb4)), "/", length(E(subb4))))
	if (l >= 5)
		print(paste("# new subgraph5 vertices / edges:", length(V(subb5)), "/", length(E(subb5))))
	if (l >= 6)
		print(paste("# new subgraph6 vertices / edges:", length(V(subb6)), "/", length(E(subb6))))
	if (l >= 3)
		print(paste("# new subgraph23 vertices / edges:", length(V(subb23)), "/", length(E(subb23))))
	if (l >= 4)
		print(paste("# new subgraph234 vertices / edges:", length(V(subb234)), "/", length(E(subb234))))
	if (l >= 5)
		print(paste("# new subgraph2345 vertices / edges:", length(V(subb2345)), "/", length(E(subb2345))))
	if (l >= 6)
		print(paste("# new subgraph23456 vertices / edges:", length(V(subb23456)), "/", length(E(subb23456))))
	
	write.graph(subb2, file= paste(name, "-subgraph-gml-2-cycles.gml", sep = ""), format = "gml")
	
	if (l >= 3)
		write.graph(subb3, file= paste(name, "-subgraph-gml-3-cycles.gml", sep = ""), format = "gml")
	if (l >= 4)
		write.graph(subb4, file= paste(name, "-subgraph-gml-4-cycles.gml", sep = ""), format = "gml")
	if (l >= 5)
		write.graph(subb5, file= paste(name, "-subgraph-gml-5-cycles.gml", sep = ""), format = "gml")
	if (l >= 6)
		write.graph(subb6, file= paste(name, "-subgraph-gml-6-cycles.gml", sep = ""), format = "gml")
	if (l >= 3)
		write.graph(subb23, file= paste(name, "-subgraph-gml-23-cycles.gml", sep = ""), format = "gml")
	if (l >= 4)
		write.graph(subb234, file= paste(name, "-subgraph-gml-234-cycles.gml", sep = ""), format = "gml")
	if (l >= 5)
		write.graph(subb2345, file= paste(name, "-subgraph-gml-2345-cycles.gml", sep = ""), format = "gml")
	if (l >= 6)
		write.graph(subb23456, file= paste(name, "-subgraph-gml-23456-cycles.gml", sep = ""), format = "gml")
}

# get cycles graph reduction
#
# cnames: cycle names
get.3cycles.subgraph.reduction <- function(cnames)
{
	g <- graph.empty(directed = FALSE)
	lc <- length(cnames)
	
	# add one vertex for each cycle
	g <- add.vertices(g, lc)
	
	# set vertices name
	g <- set.vertex.attribute(g, "name", V(g), cnames)
	
	# vname -> vids in the new graph (corresponding to cycle ids)
	vh <- hash()
	
	# vertex names (set of sorted vertex ids (e.g. "1.2.3" and "1.3.2" will both give "123"))
	svn <- c()
	
	loginfo("Computing vertices hash (vname -> vids)...")
	
	pb1 <- txtProgressBar(title = "Vertices hash computing", min = 1, max = lc, style = 3)
	
	# computing vertices hash: old vertex name -> new vids (i.e. cycle ids)
	for (i in c(1:lc))
	{
		# get ids of the vertices belonging to the cycle
		cns <- strsplit(cnames[i], "\\|")[[1]]
		
		# sort vertex ids
		# NB: works only for 3-cycles! (circular permutations)
		scns <- paste(sort(cns), collapse = "")
		
		# do not add orientations of the same set of vertices (e.g. "1.2.3" and "1.3.2")
		if (!(scns %in% svn))
		{
			svn <- append(svn, scns)
			
			# for all vertex ids of the cycle
			for (j in cns)
			{
				if (has.key(j, vh))
				{
					vids <- vh[[j]]
					vids <- append(vids, i)
					vh[j] <- vids
				}
				else
				{
					vids <- c(i)
					vh[j] <- vids
				}
			}
		}
		
		# update progress bar
		setTxtProgressBar(pb1, i)
	}
	
	close(pb1)
	
	#saveRDS(vh, "vh.rds")
	
	loginfo("Computing edges hashes (weight & vertices ids)...")
	
	pb2 <- txtProgressBar(title = "Edges hashes computing", min = 1, max = length(keys(vh)), style = 3)
	counter <- 1
	
	# ename -> weight
	ehw <- hash()
	# ename -> c(vid1, vid2)
	ehn <- hash()
	
	# draw an edge between all cycles k belongs to
	for (k in keys(vh))
	{
		vids <- vh[[k]]
		lv <- length(vids)
		
		# if the vertex k belongs to more than one cycle...
		if (lv > 1)
		{	
			for (v1 in c(1:(lv - 1)))
			{
				for (v2 in c((v1 + 1):lv))
				{
					vid1 <- vids[v1]
					vid2 <- vids[v2]
					
					ke <- paste(vid1, ".", vid2, sep = "")
					w <- ehw[[ke]]
					
					if (is.null(w))
					{
						ehw[ke] <- 1
						ehn[ke] <- c(vid1, vid2)
					}
					else
					{
						# increment the edge weight
						w <- w + 1
						ehw[ke] <- w
					}
				}
			}
		}
		
		# update progress bar
		setTxtProgressBar(pb2, counter)
		counter <- counter + 1
	}
	
	close(pb2)
	
	#saveRDS(ehw, "ehw.rds")
	#saveRDS(ehn, "ehn.rds")
	
	loginfo("Adding edges...")
	
	edges <- values(ehn, USE.NAMES=FALSE)
	weights <- values(ehw, USE.NAMES=FALSE)
	
	loginfo("Finalizing edges addition...")
	
	g <- add.edges(g, edges, weight = weights)
}

# get colinks graph
get.colink.graph <- function(core)
{
	# get cores 2 cycles
	c2 <- get.cycles(core, 2)
	
	cograph <- graph.empty(directed = FALSE)
	
	vids <- c2$v
	cids <- c2$c
	
	#print(V(core)[!(V(core)$id %in% vids)])
	
	# add vertices and labels
	cograph <- add.vertices(cograph, length(vids))
	V(cograph)$id <- vids
	V(cograph)$name <- vids
	
	indexes <- c()
	
	for (i in cids)
	{
		ncids <- as.vector(sapply(strsplit(i, "\\."), as.numeric))
		vcids <- V(core)[ncids]$id
		vindexes <- V(cograph)[V(cograph)$id %in% vcids]
		indexes <- append(indexes, vindexes)
	}
	
	# add edges
	add.edges(cograph, indexes)
}

# va: named vector of multiplicities (cycle length -> multiplicity)
# vmax: vector of the max multiplicities (cycle length -> max multiplicity)
# lv: vector of cycle lengths (>= 2)
get.score <- function(va, vmax, lv)
{
	sum <- 0
	div <- 0
	
	for (i in lv)
	{
		div <- div + exp(-i + 2)
		istr <- paste(i, sep = "")
		
		mult <- va[[istr]]
		if (is.null(mult)) next	
		sum <- sum + (mult/vmax[[istr]]) * exp(-i + 2)
	}
	
	sum/div
}

# get vertex/edge scores
#
# H: hash (vid -> named vector (cycle length -> multiplicity))
# vmax: vector of the max multiplicities (cycle length -> max multiplicity)
get.scores <- function(H, vmax, lv)
{
	scores <- c()
	
	for (k in keys(H))
	{
		score <- get.score(H[[k]], vmax, lv)
		scores[k] <- score
	}
	
	scores
}


# g: graph
# lv: vector of cycle lengths (>=2)
get.multiplicities <- function(g, lv)
{
	# hash of vertex names vectors for the different cycle lengths
	cvnames <- hash()
	# hash of edge names vectors for the different cycle lengths
	cenames <- hash()
	# hashes of edge/vertex named vector multiplicities for the different cycle lengths
	vcms <- hash()
	ecms <- hash()
	
	vmax <- c()
	emax <- c()
	
	for (i in lv)
	{
		istr <- paste(i, sep = "")
		
		loginfo(paste("Computing cycles of length: ", i, sep = ""))
		cycles <- get.cycles(g, i)
		
		cvnames[istr] <- cycles$v
		cenames[istr] <- cycles$e 
		vcms[istr] <- cycles$vcm 
		ecms[istr] <- cycles$ecm
		maxv <- max(vcms[[istr]])
		
		if (is.na(maxv))
			vmax[istr] <- 0
		else
			vmax[istr] <- maxv
		
		maxe <- max(ecms[[istr]])
		
		if (is.na(maxe))
			emax[istr] <- 0
		else
			emax[istr] <- maxe
	}
	
	# vname -> named vector (keys are in the range 2:6) and values are the corresponding vertex multiplicities
	VH <- hash()
	EH <- hash()
	
	for (i in lv)
	{
		istr <- paste(i, sep = "")
		
		# initialize named vector
		# cvnames[[istr]]: vector of vertex names belonging to cycles of length i	
		for (vid1 in cvnames[[istr]])
		{
			vidstr1 <- paste(vid1, sep = "")
			if (is.null(VH[[vidstr1]])) VH[vidstr1] <- hash()
			VH[[vidstr1]][istr] <- c()
		}
			
		for (eid1 in cenames[[istr]])
		{
			eidstr1 <- paste(eid1, sep = "")
			if (is.null(EH[[eidstr1]])) EH[eidstr1] <- hash()
			EH[[eidstr1]][istr] <- c()
		}
		
		# mult: named vector (vid -> multiplicity)
		vmult <- vcms[[istr]]
		emult <- ecms[[istr]]
		
		for (vid2 in names(vmult))
		{
			vidstr2 <- paste(vid2, sep = "")
			VH[[vidstr2]][istr] <- vmult[[vidstr2]]
		}
		
		for (eid2 in names(emult))
		{
			eidstr2 <- paste(eid2, sep = "")
			EH[[eidstr2]][istr] <- emult[[eidstr2]]
		}
	}
	
	list(vh = VH, eh = EH, vmx = vmax, emx = emax)
}

# get data frame from vertex multiplicities hash
#
# lv: set of cycle lengths
# vh: vertex multiplicities hash (vid -> (cycle length -> multiplicity))
#
# returns a data frame of vectors of lv multiplicities (one vector per vertex)
get.multiplicities.df <- function(vh, lv)
{
	df <- data.frame()
	
	# for all vertex ids
	for (i in keys(vh))
	{
		# vector of the set of lv multiplicities
		ro <- c()
		h <- vh[[i]]
		
		# for all cycle lengths
		for (j in as.character(lv))
		{
			if (j %in% keys(h))
				ro <- append(ro, h[[j]])
			else
				ro <- append(ro, 0)
		}
		
		df <- rbind(df, ro)
	}
	
	df
}

# get joint entropy of colink/clustering coeff distributions
get.graph.cc.joint.entropy <- function(g)
{
	cog <- get.colink.graph(fa)
	
	df <- data.frame()	
	
	# compute 
	for (i in V(fa))
	{
		cc <- c()
		cc <- append(cc, get.colink.coefficient(fa, cog, i))
		cc <- append(cc, get.clustering.coefficient(fa, i))
		df <- rbind(df, cc)
	}
	
	if (length(df) == 0)
	{
		0
	}
	else
	{
		# get the cross table stochastic matrix
		sm <- prop.table(table(df))
		
		get.entropy(sm)
	}
}

# get the specified graph joint entropy
get.graph.cyclic.joint.entropy <- function(g, lv)
{
	vmult <- get.multiplicities(g, lv)
	VH <- vmult$vh
	
	df <- get.multiplicities.df(VH, lv)
	
	if (length(df) == 0)
	{
		0
	}
	else
	{
		# get the cross table stochastic matrix
		sm <- prop.table(table(df))
		
		get.entropy(sm)
	}
}

##########
## testing code
##########

main <- function()
{
	# logging config
	basicConfig()
	#setLevel('ERROR')
	
	## load FA Core
	#g <- read.graph(file = "fa-core-graphml.xml", format = "graphml")
	#g <- simplify(g)
	
	## load random graph
	g <- erdos.renyi.game(20, 1/10, directed = TRUE)
	#g <- graph.ring(10)
	
	# set vertex ids
	V(g)$id <- paste(c(1:length(V(g))), sep = "")
	
	## compute cycles (Johnson)
	result <- get.johnson.cycles(g, 3)
	
	for (k in keys(result))
	{
	  print(paste("# cycles of length", k, ":", length(result[[k]])))
	  print(result[[k]])
	}
}